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SUMMARY 

The paper presents a numerical method for the solution of the conservation 
equations governing steady, reacting, turbulent viscous flow in two dimensional 
geometries, in both Cartesian and axisymmetric coordinates. These equations 
are written in Favre-averaged form and closed with a first order model. A two- 
equation K-e model, where low Reynolds number and compressibility effects are 
included, and a modified eddy-breakup model are used to simulate fluid mechan- 
ics turbulence, chemistry, and turbulence-combustion interaction. The solution 
is obtained by using a pseudo-unsteady method with improved perturbation propa- 
gation properties. The equations are discretized in space by using a finite 
volume formulation. An explicit multi-stage dissipative Runge-Kutta algorithm 
is then used to advance the flow equations in the pseudo-time. The method is 
applied to the computation of both diffusion and premix turbulent reacting 
flows. The computed temperature distributions compare favorably with experi- 
mental data. 


INTRODUCTION 

The design of future combustion concepts requires accurate numerical meth- 
ods for the simulation of liquid or gaseous fuels in a highly turbulent air- 
stream (ref. 1). These chemical reacting turbulent flows involve fluid mechan- 
ics, chemistry, and turbulence-combustion interaction. The present approach is 
devoted to the solution of the full turbulent reacting flow problem within an 
engineering accuracy, i.e., looking for a solution for averaged quantities and 
describing the reactive gas as a mixture of three species, fuel, oxidizer, and 
products, with combustion supposed to be controlled by a single step irrevers- 
ible chemical reaction. 

The reliability of a numerical method is a function of both the mathemati- 
cal modeling of the physical process and the solution algorithm. The governing 
equations are written here in Favre-averaged form. No simplification is used 
in deriving the energy conservation equation, as usually performed (ref. 8). 

The closure is achieved by using a first order model. The K-e turbulence 
model adopted here includes low Reynolds number terms, so that the equations 
are valid all over the laminar, transition, and turbulent region, as described 
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in reference 2. Furthermore, the model includes a density gradient term to 
better simulate variable density effects. The modified eddy break up mixing 
controlled reaction rate expression adopted here, introduces a concentration 
gradient dependence other than the classical concentration dependence 
(ref. 8). This model has been developed for imperfectly premixed flow 
(ref. 2), but it is useful in cases where the combustion is neither fully pre- 
mixed nor entirely diffusion controlled, as it is likely to occur under many 
ci rcumstances even if the fuel and the oxidizer enter the combustion chamber 
in separate streams. These two models allow one to simulate fluid mechanics 
turbulence, chemistry, and turbulence-combustion interaction. 

The solution is then obtained by using a pseudo-unsteady method with 
improved perturbation propagation properties (ref. 3). Pseudo-unsteady methods 
with artificial equations are quite common in turbomachinery applications 
(refs. 9 and 10), but up to now they have received little use in computing 
reacting flows, despite the fact that they are faster and simpler than classi- 
cal pressure iteration methods. At low speeds, the maximum allowable time step 
for the proposed artificial equations is indeed (ref. 12) very close to the one 
obtained in the pressure iteration method of reference 11 even by using the 
simple Lax-Wendroff algorithm. The equations are discretized in space by using 
a finite volume formulation, and integrated in the pseudo-time with an explicit 
multistage dissipative Runge-Kutta algorithm. Multi-stage schemes for the 
numerical solution of ordinary differential equations are usually designed to 
give a high order of accuracy, but in a pseudo-unsteady solution these schemes 
are selected only for their properties of stability and damping. The four 
stage adopted here allows a CFL number of about 2.6 with only negligible 
changes with reference to the simple Lax-Wendroff algorithm but with a minimum 
of artificial viscosity introduced for stability purposes. 


LIST OF SYMBOLS 

AE Arrhenius activation energy 

C constant 

c concentrations fluctuations parameter 

CFL Courant number 

D diffusion vector 

E total specific energy 

e specific internal energy 

F flux tensor 

f unknown vector 

f function 

H total enthalpy 
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Hp optimization matrix 

I identity matrix 

j index of the spatial discretization 

K turbulence kinetic energy 

k index of the multistage algorithm 

L length scale of turbulent motions 

M Mach number 

m index of the temporal discretization 

N unit outward normal 

NV updating rate 

P production of turbulence kinetic energy 

p pressure 

PF Arrhenius prefactor 

Pr Prandtl number 

Q work due to turbulence 

q heat flux vector 

R low Reynolds number term 

Re Reynolds number 

S source vector 

s stoichiometric ratio 

Sc Schmidt number 

T residual 

t time 

tu turbulence intensity 

V velocity 

v volume 

x axial coordinate 


Y mass fraction 

y radial coordinate 

(5 perturbation speed 

r diffusion coefficient 

y specific heat ratio 

AH chemical heat release 

Sr characteristic volume dimension 

St time step 

SE surface area 

e turbulence kinetic energy dissipation rate 

9 multistage scheme coefficient 

k thermal conductivity 

li viscosity coefficient 

p density 

E boundary of the fixed volume 

x viscous stress tensor 

4> reaction rate 

4) Schwab-Zel dovi ch function 

Q artificial viscosity coefficient 

Subscripts 
i inviscid 

fu fuel 

1 laminar 

ox oxidizer 

t turbulent 

v viscous 
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Governing Equations 


The unknown vector f is the solution of a system of conservation equa- 
tions. This system is written in Favre-averaged, dimensionless, vector, inte- 
gral form as fol lows 


Ij. J N ■ (F I * F v>« - J J v J 5 d * 

The basic dependent variables are density, velocity and energy, 
conservation equations read as follows 


Their 


fp • V 

F. - p • V V + p • I 

\p * E • V + p • I • Vy 


f » ■ £ . V ♦ J 

(0 
s = I o 

Vp • e - P + 0 + AH • <j>^ u . 


where from the equation of state of a perfect gas 

p = ( Y - 1 ) • p • e 


wi th 


e = (E - V 2 / 2) = H - V 2 /2)/y 

The reactive variables are the fuel mass fraction and a Schwab-Zeldovi ch 
function 

* * Y fu - V s 

Their conservation equations are written as follows 
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The stress tensor x is the sum of a laminar and a turbulent part, where 
the latter is expressed according to a classical eddy viscosity concept, i.e., 

t = - | • (ji/Re • div(V) + p • K) • I + 2 . p/Re • (def(V)) 


where 


Similarly, the heat flux vector is given by 

q = -y • k • grad(e)/(Re • Pr) 
while the diffusion vector is given by 

Dy = -Ty • grad(Y)/(Re • Sc y ) 

where 


K = 


r - r a + r t 


The turbulent viscosity coefficient is expressed accordingly with the 
Prandtl -Kolmogorov formulation 




= C • f • p • Re 
P P 


K 2 /e 


while 


T ~ ~ p 

The turbulence variables are the turbulence kinetic energy and its dissi- 
pation rate. Their conservation equations are written in the following low 
Reynolds number and compressible form 



i 
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s = 


'P - 0 - p 


+ R. 


(P - 0) • C J • f„, • e/K - C 


'e2 


f 


e2 


e/K + R 


where the production of turbulent energy from the mean flow energy P 
work due to turbulence Q are given by 


and the 


* 


(def(V) : def(V) )/Re 


(grad(p) -grad(p>/p 2 )/Re 


v ' v p H t 

the low Reynolds number functions are given by reference 2 

f = exp(-3.4/( 1 + Re t /50)) 


f e] - i.o 


f 2 =1.0- 0.33*exp(-Re 3 ) 


R = -2 
e 




(grad(e) ) /Re 


R k .-2 


(gradOC) ) <: /Re 


where 


Re. = p 


Re/(p c 


e ) 


and the model constant are assumed as 


c e1 “ ] ' 43 


C e2 = 1.92 


Sc = 1.3 

e 


Sc K = 1 .0 


C = 0.09 


C = 1.0 
p 


The above two equation models do not take into account the preferential 
damping of velocity fluctuations in the direction normal to the wall, but it is 
quite general and it is useful in laminar, transition and turbulent regions. 
Furthermore, the model adopts the assumptions and approximations which are nor- 
mally used for constant density flows, by retaining the gradient diffusion 
model to be rewritten in the density weighted form without any explicit account 
being taken of density fluctuations. However, the introduction of the compres- 
sibility term Q allows a partial consideration of variable density effects. 

The reaction rate is finally expressed according to a modified eddy- 
breakup model. The mixing controlled reaction rate is given as follows 


^fu.mix "Si p 


. 1/2 


K /( 
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where for imperfect premixing we may suppose 

c = c(Y fu> <j>, grad(Y fu >, grad(<t>)) 

We use here (ref. 2) 

c - min(y2 u , (Y ox /s) 2 , K 3 /, 2 • grad<Y fu > 2 • C +2 ) 

The kinetics controlled reaction rate is then evaluated as follows 

Vkin " p2 ' V fu ' Y ox • PF ' 

Finally, the actual reaction rate is taken to be the smaller of the values 
from the previous expressions. 

Along the inflow boundary, all the flow variables except the static pres- 
sure are specified. Along the outflow boundary, only the static pressure is 
specified. Along the solid boundaries, the no-slip condition requires the van- 
ishing of velocity, turbulence kinetic energy, and turbulence kinetic energy 
dissipation rate, the latter being obviously intended as the modified quantity 
used in the conservation equations. The specific internal energy is then set 
equal to the value following from the wall temperature for constant tempera- 
ture walls, or the energy gradient normal to the wall is set equal to zero for 
adiabatic walls. Furthermore, the species concentration gradient normal to 
the wall is set equal to zero. Along a symmetry plane (axis) the normal 
derivatives of all the flow parameters vanish, with the exception of the normal 
velocity component, to be set equal to zero. 


NUMERICAL SOLUTION 

The proposed equations are solved in two-dimensional geometries. The der- 
ivation of the two-dimensional equations from the previous general form is 
straightforward and not presented here. The equations are written in a form 
useful in both Cartesian and axisymmetric geometries by introducing a switch 
parameter distinguishing between these two coordinates (ref. 2). 

The discretization adopted is a pseudo unsteady, finite volume, dissipa- 
tive, explicit one. In the proposed pseudo unsteady, the solution of the 
steady equations is obtained as the asymptotic solution of the following arti- 
ficial unsteady equations 



■* 

m . n-'dv ♦ 

■> 

' N • (F. + F )dE = 

> 

-> 

J 

V J 

V3ty p J 

E J 

1 v J 

U 

V J 


These unsteady equations are generally constructed in order to obtain the bet- 
ter convergence rate, obviously providing that the steady state solution is 
not al tered. 

From the identity between the convergence process and the elimination 
process of the initial perturbations to the steady solution, the convergence 


i 
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parameters are determined in order to improve the perturbation propagation or 
damping. We use (ref. 3) 


0 . . . 0 ; 

0 . . . 0 

0 . . . 0 

f 2 . . . 0 

0 ... 1 

f 1 = min(M^ -1,0) 
f 2 = max(f 1 + 1 , C > 

wi th 

M = V . < Y • p/p ) -1 /2 
Cp is a small positive number, and 

B max ' min<l3 P B 2> 

8, - V • {<1 ♦ M 2 )/2 + [(() - M 2 >/2> 2 ♦ 1] ,/2 } 


where 


H P “ 6r/( W 


1 

0 

0 

V 2 / 2 


0 

1 

0 


0 

0 

1 


- f 1 


V x - *1 


P 2 - V • <1 + 1/M) 

This produces an improved propagation, but only for subsonic flows. 

These equations are then discretized in space by using a finite volume 
discretization. The mesh is nonorthogonal and curvilinear, conforming to the 
boundaries of the domain, with lines intersecting at arbitrary angles, properly 
refined where high gradients are expected to occur. The discretization nodes, 
located at the intersection of these lines, are the centers of hexagonal con- 
trol volumes, obtained by connecting the six surrounding nodes. A sample com- 
putational domain and the hexagonal control volume are shown in figure 1. 

The discretized equations are written as follow (ref. 2 and 4) 


3 £ _ 

8t " 



' 6 
J-l 


(F i.j + F v,j ),N j' SE j 


v + S 


where the subscript j refers to every face of the finite volume. 
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The equations are finally discretized in time by using an explicit, dis- 
sipative discretization. Let the previous equation be rewritten with the addi- 
tion of a dissipative term as follows 

|| = T(f ) + D(f ) 

where T represents the residual and D is the dissipative term. The 
explicit k-stage Runge-Kutta algorithm is written as follows (St = CFL) 

^(0) _ ^m 



The subscript j refers now to every surrounding node involved in the 
definite volume approximation, and the superscript m refers to local time 
nvSt. The terms referring to time m*-St are updated only at specific itera- 
tions and assumed constant between two updatings. The updating rate is taken 
equal to NV iterations, to be determined accordingly with a numerical 
optimization. 

A four-stage scheme, with the standard coefficients 

®] - 5 e 2 * 5 e 3 * 1 

has a Courant limit of about CFL = 2.6. 

The time step is evaluated according to the classical CFL stability 
limit all over the computational domain. It is taken slightly smaller than the 
local CFL number in order to take into account the neglected stability limit 
due to the viscous terms. 
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Thanks to the efficient pseudo unsteady solution, the method appears to be 
rather fast, while the explicit finite volume discretization allows ease of 
understanding and computer programming. 


RESULTS 

First calculations were performed for a partially premixed turbulent 
reacting flow. Figure 2 shows schematically the flow domain. Experimental and 
theoretical results are presented in reference 13. In this configuration, mix- 
ing of two parallel streams, one of hot gases and the other of a fresh mixture 
of air and methane, in a constant area duct is considered. The hot jet causes 
a flame to be ignited and stabilized in the duct. The inlet duct, with a cross 
section of 100 by 100 mm 2 , is split into two parts. The upper section (80 by 
100 mm 2 ) is assigned to the fresh air and methane mixture, with a velocity of 
65 m/s, a temperature of 580 K, and a mixture ratio of 0.8; the lower one 
(20 by 100 mm 2 ) is assigned to the pilot flame made up of hot gases, with a 
velocity of 130 m/s and a temperature of 2000 K; the walls are insulated. 

Calculations were performed on a 35 by 45 computational grid, nonuniform 
but orthogonal. Conventional steady solutions were obtained in about 2x1Ch 
iterations with a total CPU time of about 1-1/2 hr on a VAX 8800. Figures 3 
and 4 show a comparison of predicted /method 2/ and measured transverse temper- 
ature profiles at x = 42 mm and at x = 122 mm. The agreement in both sec- 
tions lies within engineering accuracy. The introduction of the influence of 
kinetics in the evaluation of the reaction rate and a better calibration of the 
model constants leads to a better accuracy than the one obtained in reference 2 
/method 1/, even if at x = 122 mm the temperatures in the mixing layer are 
still underestimated. 

The method has then been applied to a turbulent diffusion reacting flow. 
Figure 5 shows schematically the flow domain. Experimental and theoretical 
results are presented in reference 8. In this configuration, there is a cen- 
tral jet of fuel, substantially methane, with a velocity of 21.3 m/s and a tem- 
perature of 300 K, and a concentric jet of oxidizer, with a velocity of 
34.3 m/s and a temperature of 589 K. The outer radius of the fuel jet is 8 mm, 
the inner and outer radius of the air jet are respectively 11.1 and 28.6 mm. 

The radius of the combustor is 101.6 mm. The wall temperature is 1140 K. 

Calculations were performed on a 45 by 45 computational grid, nonuniform 
but orthogonal. Conventional steady solutions were obtained in less than 3x10^ 
iterations with a total CPU time of about 3 hr on a VAX 8800. Figures 6 and 7 
show a comparison of predicted /method 2/ and measured radial temperature prof- 
iles at x = 95 mm and x = 398 mm. 

Even if the accuracy appears to be not optimal, the proposed method allows 
a better accuracy than the method presented in reference 8 /method 1/, despite 
the single step chemical reaction adopted here. This is substantially due to a 
better modeling of the reaction rate. The mixing controlled reaction rate of 
reference 8 is simply taken proportional to the smaller value of the fuel and 
the oxidizer concentration. The present mixing controlled reaction rate is 
taken proportional to the fuel concentration gradients with limitations arising 
from the availability of fuel and oxidizer. 
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In a laminar diffusion flame, little oxidizer is detectable within the 
reaction zone envelope, and little fuel is detectable outside of it. The reac- 
tion zone is a very thin envelope, and fuel and oxidizer concentrations show a 
small overlapping. In a turbulent diffusion flame, the reaction zone is 
thicker, and the averaged values of fuel and oxidizer overlap significantly. 

The reaction zone is now a finite volume, even if always relatively small. The 
introduction of the gradient concentration dependence other than the concentra- 
tion dependence appears to lead to a better simulation of the thickness of the 
reaction zone. 


CONCLUSIONS 

The paper has presented a numerical method for the study of steady, react- 
ing turbulent viscous flow in two-dimensional geometries, both Cartesian and 
axi symmetric . 

The proposed two-equation K-e model, where low Reynolds number and com- 
pressibility effects are included, and a modified eddy-breakup model give a 
satisfactory description of fluid mechanics turbulence, chemistry, and 
turbulence-combustion interaction within an engineering accuracy. 

The pseudo-unsteady method with improved perturbation propagation proper- 
ties and the explicit multi-stage dissipative Runge-Kutta algorithm appear to 
be fast and rel iable. 

The application of the method to the computation of the temperature dis- 
tributions for a diffusion and a premixed turbulent reacting flow shows a 
satisfactory agreement between experimental and computational results. 
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FIGURE 1. - SAMPLE COMPUTATIONAL DOMAIN AND 
HEXAGONAL CONTROL VOLUME. 
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FIGURE 2. - PREMIXED TURBULENT REACTING FLOW: PHYSICAL DOMAIN. 
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FIGURE 3. - PREMIXED TURBULENT REACTING FLOW: COMPARISON 
OF PREDICTED AND MEASURED TEMPERATURE DISTRIBUTIONS. 
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FIGURE 4. - PREMIXED TURBULENT REACTING FLOW: COMPARISON 
OF PREDICTED AND MEASURED TEMPERATURE DISTRIBUTIONS. 
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FIGURE 5. - DIFFUSION TURBULENT REACTING FLOW: PHYSICAL DOMAIN. 
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FIGURE 6. - DIFFUSION TURBULENT REACTING FLOW: COMPARISON 
OF PREDICTED AND MEASURED TEMPERATURE DISTRIBUTION. 
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FIGURE 7. - DIFFUSION TURBULENT REACTING FLOW: COMPARISON 
OF PREDICTED AND MEASURED TEMPERATURE DISTRIBUTIONS. 
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